﻿#include <occmesh/api/meshfit.h>

MeshFit::MeshFit(MeshPtr mesh) : m_mesh(mesh)
{
}

TopoDS_Face MeshFit::Build(std::vector<std::vector<int>> &cornerNodes, std::vector<std::vector<int>> &boundaryNodes)
{
    MeshPtr mesh = m_mesh;
    std::cout << "Selected mesh..." << std::endl;
    std::cout << "Vertices: " << mesh->n_vertices() << std::endl;

    Parameterize();

    auto surface = BuildSurface();
    if (cornerNodes.empty() || cornerNodes[0].size() < 3)
        return TopoDS_Face{};

    FixNodes(cornerNodes[0], boundaryNodes[0]);
    auto curves = BuildBoundaryCurves(cornerNodes[0], boundaryNodes[0]);
    BRepBuilderAPI_MakeFace faceBuilder(surface, OCCUtils::TrimmSurface(surface, curves));

    for (int i = 1; i < cornerNodes.size(); ++i)
    {
        if (cornerNodes[i].size() < 3)
            continue;

        for (int i = 0; i < cornerNodes.size(); ++i)
            FixNodes(cornerNodes[i], boundaryNodes[i]);

        auto curves = BuildBoundaryCurves(cornerNodes[i], boundaryNodes[i]);
        faceBuilder.Add(OCCUtils::TrimmSurface(surface, curves));
    }

    std::cout << "Finish" << std::endl;
    return faceBuilder.Face();
}

void MeshFit::Parameterize()
{
    MeshPtr mesh = m_mesh;

    TIMER_START(LSCM);
    // SODFeatureMetric metric(15);
    // LSCMBuilder lscmBuilder(mesh.get(), &metric);
    LSCMBuilder lscmBuilder(mesh.get());
    auto charts = lscmBuilder.Build();
    TIMER_END(LSCM);

    TIMER_START(ARAP);
    ARAPBuilder arapBuilder(mesh.get(), charts, 1e-16);
    arapBuilder.Build(3);
    TIMER_END(ARAP);
}

Handle(Geom_BSplineSurface) MeshFit::BuildSurface()
{
    MeshPtr mesh = m_mesh;
    TIMER_START(BuildSurface);

    std::vector<gp_Pnt> cloud;
    std::vector<gp_Pnt2d> uv;
    float umax = -std::numeric_limits<float>::max();
    float vmax = -std::numeric_limits<float>::max();
    float umin = std::numeric_limits<float>::max();
    float vmin = std::numeric_limits<float>::max();
    for (auto vh : mesh->vertices())
    {
        auto point = mesh->point(vh);
        auto coord = mesh->texcoord2D(vh);
        cloud.push_back(gp_Pnt(point[0], point[1], point[2]));
        uv.emplace_back(coord[0], coord[1]);
        umax = std::max(umax, coord[0]);
        vmax = std::max(vmax, coord[1]);
        umin = std::min(umin, coord[0]);
        vmin = std::min(vmin, coord[1]);
    }
    double ratio = (vmax - vmin) / (umax - umin);
    std::cout << "Ratio: " << ratio << std::endl;

    int Nu = 20, Nv = Nu * ratio;
    if (Nv < 4)
    {
        Nv = 4;
        Nu = Nv / ratio;
    }
    int p = 3, q = 3;
    std::vector<double> w(cloud.size(), 1.0); // 每个点的拟合权重
    double alpha = 1e-6;                      // 弦长最小化权重
    double beta = 1e-3;                       // 二阶平滑化权重
    auto surface = OCCUtils::Chord(cloud, uv, Nu, Nv, p, q, w, alpha, beta);

    TIMER_END(BuildSurface);

    return surface;
}

void MeshFit::FixNodes(std::vector<int> &cornerNodes, std::vector<int> &boundaryNodes)
{
    TIMER_START(FixNodes);

    // 注意 boundaryNodes 的顺序可能与 cornerNodes 相反，需要首先检测一下
    int recordPos[] = {0, 0};
    for (int i = 0; i < boundaryNodes.size(); ++i)
    {
        if (boundaryNodes[i] == cornerNodes[1])
            recordPos[0] = i;
        else if (boundaryNodes[i] == cornerNodes[2])
            recordPos[1] = i;
    }
    if (recordPos[0] > recordPos[1])
    {
        cornerNodes.push_back(cornerNodes[0]);
        std::reverse(cornerNodes.begin(), cornerNodes.end());
        cornerNodes.pop_back();
    }

    TIMER_END(FixNodes);
}

std::vector<Handle(Geom2d_Curve)> MeshFit::BuildBoundaryCurves(const std::vector<int> &cornerNodes,
                                                               const std::vector<int> &boundaryNodes)
{
    MeshPtr mesh = m_mesh;
    TIMER_START(Trimm);

    int nextCornerId = 1;
    int prevBoundaryId = 0;

    OpenMesh::VertexHandle startvh(cornerNodes[0]);
    std::vector<Handle(Geom2d_Curve)> curves;
    for (int i = 0; i < boundaryNodes.size(); ++i)
    {
        if (boundaryNodes[i] != cornerNodes[nextCornerId])
            continue;

        std::vector<gp_Pnt2d> tuv;
        for (int j = prevBoundaryId; j <= i; ++j)
        {
            OpenMesh::VertexHandle vh(boundaryNodes[j]);
            auto coord = mesh->texcoord2D(vh);
            tuv.push_back(gp_Pnt2d(coord[0], coord[1]));
        }
        auto curve2d = OCCUtils::Interpolate(tuv);
        curves.push_back(curve2d);
        nextCornerId = (nextCornerId + 1) % cornerNodes.size();
        prevBoundaryId = i;
    }

    gp_Pnt2d corners[3];
    corners[0] = curves[0]->Value(curves[0]->FirstParameter());
    corners[1] = curves[1]->Value(curves[1]->FirstParameter());
    corners[2] = curves[2]->Value(curves[2]->FirstParameter());

    // 检查方向并修正
    // bool isAntiClockwise = corners[0].X() * (corners[1].Y() - corners[2].Y()) +
    //                        corners[1].X() * (corners[2].Y() - corners[0].Y()) +
    //                        corners[2].X() * (corners[0].Y() - corners[1].Y()) >
    //                    0;
    // if ((antiClockwise && !isAntiClockwise) || (!antiClockwise && isAntiClockwise))
    // {
    //     std::reverse(curves.begin(), curves.end());
    //     std::transform(curves.begin(), curves.end(), curves.begin(), [](Handle(Geom2d_Curve) curve) {
    //         return curve->Reversed();
    //     });
    // }

    TIMER_END(Trimm);

    return curves;
}